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Abstract 

We compute the electric dipole moment d n of the neutron from a fully dynamical simula¬ 
tion of lattice QCD with 2+1 flavors of clover fermions and nonvanishing 9 term. The latter 
is rotated into a pseudoscalar density in the fermionic action using the axial anomaly. To 
make the action real, the vacuum angle 6 is taken to be purely imaginary. The physical value 
of d n is obtained by analytic continuation. We find d n = -3.9(2)(9) x 10 -16 6 e cm, which, 
when combined with the experimental limit on d n , leads to the upper bound \9\ < 1 Ax I (T 11 . 
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1 Introduction 


The electric dipole moment d n of the neutron provides a unique and sensitive probe to physics 
beyond the Standard Model. It has played an important part over many decades in shaping and 
constraining numerous models of CP violation. While the CP violation observed in K and B 
meson decays can be accounted for by the phase of the CKM matrix, the baryon asymmetry of 
the universe cannot be described by this phase alone, suggesting that there are additional sources 
of CP violation awaiting discovery. 

QCD allows for CP-violating effects that propagate into the hadronic sector via the so-called 
6 term S e in the action, 

S = S o + S o, S g = i 6 Q , (1) 

where (in lattice notation) 

Q = a' £ f;,f € z ( 2 ) 

is the topological charge, and S 0 is the standard CP-preserving QCD action. Thus, there is 
the possibility of strong CP violation arising from a nonvanishing vacuum angle 6. In a wide 
class of GUTs the diagrams that generate a high baryon to photon asymmetry contribute to 
the renormalization of 6, and hence to the electric dipole moment of the neutron. With the 
increasingly precise experimental efforts to observe the electric dipole moment [[H (21 HI, it is 
important to have a rigorous calculation directly from QCD. 

It is practically impossible to perform Monte Carlo simulations with the action ([[]) in four 
dimensions for any sensible definition of the topological charge and any angle |0| > 0. Absorbing 
the 6 term into the observable |4l[5]| is not a viable alternative, as (Q 2 ) is found not to vanish if 
one of the quark masses is taken to zero at present values of the coupling. In Fig. Q] we show the 
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Figure 1: The topological susceptibility on the SU(3) symmetric line m u = m d = m s as a function 
of ml in units of to. 
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topological susceptibility^ = ( Q 2 )/V on 32 3 x64 lattices taken from @ at spacing a = 0.074 fm. 
The charge Q has been computed from the Wilson flow 0 at flow time to- Similar results have 
been reported in 0. As a result, d n will not vanish in the limit of zero quark mass either, 
except perhaps for chiral fermions. Exactly that was found in 0. This precludes a meaningful 
extrapolation of d n to the physical point. There are indications that the situation will improve for 
lattice spacings a < 0.04 fm only [|E]- 

It so happens that the 9 term can be chirally rotated into the fermionic part of the action, 
making use of the axial anomaly IflOll . The outcome of that is 

S e = - l -6 m a A ^ ( uy 5 u + dy 5 d + vy 5 ,y) , m~ l = ^ ( m~ l + m~ d l + m~'^j (3) 

* 

for three quark flavors with nondegenerate masses. This action lends itself to numerical simula¬ 
tions for imaginary values of 9 IfTTTl . As we are mainly interested in small values of 9, the results 
can be analytically continued to real numbers without difficulties, assuming that the theory is 
analytic in the vicinity of 9 = 0. 

In this paper we present an entirely dynamical calculation of the electric dipole moment of 
the neutron on the lattice. This is a challenging task. As d n quickly diminishes towards physical 
quark masses, the angle 9 has to be chosen increasingly larger to compensate for that. This in turn 
leads to a substantial increase of zero modes, which slows down the simulations substantially and 
eventually will result in exceptional configurations Ill2l . 


2 The simulation 


We follow mm an d start from the SU(3) flavor symmetric point m u = m d = m s = m 0 , where 
m„ = m K . Our strategy has been to keep the singlet quark mass m = (m u + m d + m s )/3 fixed 
at its physical value, while 6m q = m q - m is varied. As we move from the symmetric point to 
the physical point along the path m = constant, the 5 quark becomes heavier, while the u and d 
quarks become lighter. These two effects tend to cancel in any flavor singlet quantity, such as the 
topological susceptibility^, = (Q 2 )/V. The cancellation is perfect at the symmetric point [j6]|. 

We assume u and d quarks to be mass degenerate, writing m ( = m u = m d . The vacuum angle 
is taken purely imaginary, 

9 = id. (4) 

This leads us to consider the action 


Se = 9 


m, m s 
2 m s + mt 


y. (wy 5 u + dy 5 d + sy 5 s} 


(5) 


which is real and vanishes at m c = 0 as well as m. v = 0. 

Our fermion action has single level stout smearing for the hopping terms together with un¬ 
smeared links for the clover term. With the (tree level) Symanzik improved gluon action this 
constitutes the Stout Link Non-perturbative Clover or SLiNC action [|T4l . To cancel O(a) terms 
the clover coefficient c S w has been computed nonperturbatively. For each flavor the fermion 
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# 

K ( 

Ks 

am n 

am K 

am N 

A 

1 

0.12090 

0.12090 

0.1747(5) 

0.1747(5) 

0.4673(27) 

0.003 

2 

0.12090 

0.12090 

0.1747(5) 

0.1747(5) 

0.4673(27) 

0.005 

3 

0.12104 

0.12062 

0.1349(5) 

0.1897(4) 

0.4267(50) 

0.003 

4 

0.12104 

0.12062 

0.1349(5) 

0.1897(4) 

0.4267(50) 

0.005 


Table 1: The simulation parameters with m = constant. The hadron masses refer to A = 0. 


action to be simulated reads 


S* = ^ + ^ = a 4 

X ' 


1 

w [av Ffiv 


+ + Ta ys ) 


q. 


( 6 ) 


where D is the Wilson Dirac operator and 


A = 62a 


ni/ m s 

2 m s + in/ 


(7) 


The extra term in the action © can be treated in a similar way as we treat disconnected diagrams 
in calculations of singlet hadron matrix elements and renormalization factors llT5l IT6ll . We use 
BQCD lfT7ll to update the gauge fields. The calculations are done on 24 3 x48 lattices at /3 = 5.50. 
At this coupling the lattice spacing was found to be a = 0.074(2) fm lfT8l . using the center of 
mass of the nucleon octet to set the scale. The parameters of the simulations are listed in Tabled] 
Each ensemble consists of (7(2000) trajectories. The quark masses on the m = constant line are 
given by m q = l/2/q y - 1 /2 a- 0 ,c- with a f = 0.12110 |f6l. 

We expect our ensembles to carry nonvanishing topological charge, ( Q) oc -6(Q 2 ) C , with 
(Q 2 ) c = (Q 2 ) ~ (Q) 2 00 in ITU?! . In Fig. [2] we show the charge histogram for ensemble 4, together 
with a Gaussian fit. As before, the topological charge has been computed from the Wilson flow 
at flow time to [H. Evidently, Q peaks at negative values. In Fig. [3] we show ( Q) as a function 
of 6 for both sets of quark masses, together with linear plus cubic fits. We find the slopes of the 
individual curves to be approximately proportional to m, as expected. 


3 The evaluation 


At nonvanishing vacuum angle 6 the nucleon matrix element of the electromagnetic current reads 
in Euclidean space 

(p\ s'\J/i\p, s) = u g (p', s')Jg ue(p, s ), ( 8 ) 

where 


g 2 E"(r) 2 g 2 F^r) 

J/1 = VgFM ) + CTfjvC/v' ~ — + (yq q tJ - qr)y 5 F A (q ) + cr MV q v y 5 — 


2m! 


(9) 


N N 

and q - p' - p, q 2 = (p' - p) 2 - ( E dr - E 9 ) 2 . In the 6 vacuum the Dirac spinors pick up a 
phase 


Ug(p, s ) = e ,aW)y5 u(p, s ), 
Ug(p, s ) = u(p, s ) e ,a(9)y5 , 
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Figure 2: The topological charge distribution p{Q) (with p(Q ) = 1) of ensemble 4 at Kc = 

0.12104, k s = 0.12062 and A = 0.005, together with a Gaussian fit. 



Figure 3: The average charge (Q) as a function of 6 for ensembles 1 and 2 (#) and ensembles 3 
and 4 (■), together with linear plus cubic fits. 
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Figure 4: Disconnected insertion of the pseudoscalar density to lowest order. Gluon lines are 
omitted. 


so that 

2 uo(p, S) Up, s) = Q ia(ff)ys | -^+ < j eM eyy 5 ( 11 ) 

with yp = yp + iEy 4 . The electric dipole moment is given by 


eF°( 0 ) 

2m N 


( 12 ) 


The topological 6 term ([[]) polarizes the vacuum. Diagrammatically it solely contributes to 
internal gluon lines. Similarly, the flavor-singlet pseudoscalar density in ([5]) and © interacts with 
the nucleon through quark-line disconnected diagrams only [STL 221- This is sketched in Fig.Hl 
Consequently, the quark propagators in the nucleon matrix element ® are computed with the 
action S',, neglecting the S q 0 term. 

We denote the two-point function of a nucleon of momentum p in the 6 vacuum by G e NN (t, p). 
The phase factor a is obtained from the ratio of two-point functions 


Tr [G e NN (t; 0)C 4 ] 
Tr [G e NN (t\ 0)C 4 y 5 ] 


1+C< f ttOT ;l Z/e-*. 


(13) 


where r 4 = (1 + y 4 )/2. Equation ([T3T) defines m e N , the nucleon mass for the action ©, and Z N . 
The form factor F 3 (q 2 ) can be extracted from the ratio of three-point and two-point functions, 
generalizing the methods developed in Il23l 


Rn(t', t; p', p) 


G n\n^^P',P) 
Tr \G° NN (t'\p')F 4 ] 


Tr [G e NN (f, p')r 4 \ Tr [G e NN (t'; p')r 4 ] Tr [G 0 NN (f - f; p)r 4 \ \ 1/2 
Tr [G e NN {t- p)r 4 \ Tr [G e NN (f, p)r 4 \ Tr [G e NN (f - t\ p')r 4 \ J 


( 14 ) 


E 6, E B 


(E f >' + ml) (E B + m° N ) 
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Figure 5: The phase factor a(6) as a function of 6 for our two sets of quark masses. 


where G° h [j N (t', t; p',p ) is the three-point function, with t' being the time location of the nucleon 
sink and t the time location of the current insertion, and the function F(T, J~ M ) is 


f(U„) = 7 Trr 


. E e 'y 4 - iyp' + m e N . 

ia( 0)y 5 - ' > _ N_ e ia(6)y 5 

E er 




E 6 


(15) 


with given in ([9]). The three-point functions are calculated for various choices of nucleon 
polarization, / = f 4 , ir 4 y 5 y u iE 4 y 5 y 2 and ir 4 y 5 y 3 . For J fl we take the local vector current qy^q. 


4 Results 

In physical units, the pion and kaon masses areQ 


Kt 

K.s 

m n [MeV] 

m K [MeV] 

0.12090 

0.12090 

465(13) 

465(13) 

0.12104 

0.12062 

360(10) 

505(14) 


(16) 


To a good approximation 2 mr K + ml = constant, in accord with the leading order chiral expansion 
2 m 2 K + ml = 6 Bq in. 

At imaginary values of 6, both a{0) and F° 3 are imaginary. Thus, we can write 

a(0) = i a(6 ), F\ = i F\ . (17) 

In Fig.[5]we show the results for the phase factor of 60, and in Fig. [6] we show the form factor Fy" 

'it is to be noted that the pseudoscalar mass at our flavor symmetric point are somewhat larger than the physical 
value + m 2 K+ + m /3 =413 MeV. 
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Figure 6: The ratio of form factors F°’ n /F S l ’ p for K/ = k s = 0.12090 and A = 0.005. 



Figure 7: The renormalized form factor F R,nR (0) as a function of 6, together with a linear plus 
cubic extrapolation, Ff n R (0 ) = A 6 + B 6 3 , to 6 = 0. 
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Figure 8 : The dipole moment of the neutron extrapolated to the physical point along the path 
m = constant. 

of the neutron divided by Fof the proton for ensemble 2. If the radii of the two form factors 
are close to one another, the q 2 dependence is largely cancelled out in the ratio. Indeed, the 
ratio shows only a mild q 2 dependence and thus may be extrapolated linearly to q 2 = 0. The 
extrapolated value is the renormalized form factor Ff" R (0), using the fact that F A ' pR (0) = 1, 
from which we obtain the electric dipole moment (fl2l) . In Fig.[7]we show our results for F'l" A ’(0) 
as a function of 0 for our two sets of quark masses. It should be noted that the actual expansion 
parameter is A, given in ©, which is a very small number. 

Ultimately, we are only interested in F®(0) (we drop the superscripts n, R on F 3 from now on) 
at very small values of 0. Even so, we do not have sufficient data to constrain the extrapolation 
of F®(0) to 0 = 0. This will result in a systematic error. To estimate the error, we have employed 
a linear plus cubic fit, A 0 + B 0 3 , a Pade fit, A 0/(1 + B 0 2 ), allowing for corrections of O(0 5 ) and 
higher, as well as a linear fit, A 0, to the lowest 6 point each. We identify the central value of A 
with the derivative of F|(0) at 0 = 0, F^ (0). The coefficient A of the linear plus cubic fit shown 
in Fig. [7] turns out to be close to the central value. The error of F^\ 0) is estimated to be the 
largest deviation of A from the central value. After continuing 0 and F®(0) back to real values, 
we finally obtain, writing d n = e F? ’(()) 6/2m N , 


m n [MeV] 

m K [MeV] 

d n [e fm 0] 

465(13) 

360(10) 

465(13) 

505(14) 

-0.0297(38) 

-0.0215(25) 


To extrapolate (fl8l) to the physical point, we make use of the analytic expressions derived from 
covariant U{3) L x U(3) R baryon chiral perturbation theory in Il24l to NLO, with the additional 
constraint 2 m 2 K + m 2 = constant oc m. This basically involves one free low-energy constant, 
w a (jj), only. A fit to the lattice data gives w a (/u = 1 GeV) = 0.04(1) GeV~'. The result of the fit 
is shown in Fig. [U Note that d n vanishes at 2 m 2 K - ml = 0 due to the constraint m = constant. At 
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the physical point this finally leads to 


d n = -0.0039(2)(9) [e fin 6 ]. (19) 

The first error is purely statistical. The second error is a conservative estimate of NNLO effects. 
It covers the naive result from a polynomial extrapolation, d n = -0.0043 [e I'm ()\. 

Our result (fl9l) translates into constraints on CP violating contributions to the action at the 
quark and gluon level. The current experimental bound on the electric dipole moment of the 
neutron is ll25ll \d n N \ < 2.9 x 10 -13 [e fm]. Combining this bound with (fl9l) . we arrive at the upper 
bound on 6, 

\e\ < 1A x 10 -11 . (20) 


5 Conclusions 

It should be noted that in this exploratory work we have not included contributions from dis¬ 
connected insertions of the electromagnetic current. However, since these contributions vanish 
exactly at the flavor symmetric point, we do not expect them to have a significant effect to our 
conclusions. It remains to be seen how big they are at the physical point. 

The vacuum angle 6 renormalizes as 0 s = (Z| /Z P ) 6, where Z s s and Z P are the renormaliza¬ 
tion constants of the flavor-singlet scalar density and the pseudoscalar density, respectively. In 
the continuum Z s s /Z P = 1. A caveat of our calculations is that clover fermions, though 0(a) 
improved, break chiral symmetry at finite lattice spacings. On our present lattices Z s s /Z P = 
0.8 - 0.9 ll26l l6l[T6l. which might imply a systematic error of 0(10%). 

To sum up, we have successfully computed the electric dipole moment of the neutron d n from 
simulations of 2 + 1 flavor lattice QCD at imaginary vacuum angle 6, using the axial anomaly to 
rotate the topological charge density into a flavor singlet pseudoscalar density in the fermionic 
action. Only disconnected insertions of the pseudoscalar density contribute to the dipole moment, 
which required the generation of new gauge field ensembles with the modified action ©. Clearly, 
our results will have to be substantiated by simulations on larger lattices, at smaller pion masses 
and smaller lattice spacings, as well as for a wider range of A parameters. This is a challenging 
task, which we hope to report on in due course. 
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